#include <math.h>
#include <vector>
#include <list>
#include <stdio.h>
#include <algorithm>

using namespace std;

template<int dim> struct Point
{
    double x[dim];
    Point(const Point& p) {
        for (int i = 0; i < dim; i++) x[i] = p.x[i];
    }

    Point& operator= (const Point &p) {
        for (int i = 0; i < dim; i++) x[i] = p.x[i];
        return *this;
    }
    bool operator== (const Point &p) {
        for (int i = 0; i < dim; i++) if (x[i] != p.x[i]) return false;
        return true;
    }
    Point(double x0 = 0.0, double x1 = 0.0, double x2 = 0.0) {
        x[0] = x0;
        if (dim > 1) x[1] = x1;
        if (dim > 2) x[2] = x2;
        if (dim > 3) throw("Point not implemented for DIM > 3");
    }
};

struct Circle {
    Point<2> center;
    double radius;
    Circle(const Point<2> &cen, double rad) : center(cen), radius(rad) {}
    };
    Circle circumcircle(Point<2> a, Point<2> b, Point<2> c) {
    double a0,a1,c0,c1,det,asq,csq,ctr0,ctr1,rad2;
    a0 = a.x[0] - b.x[0]; a1 = a.x[1] - b.x[1];
    c0 = c.x[0] - b.x[0]; c1 = c.x[1] - b.x[1];
    det = a0*c1 - c0*a1;
    if (det == 0.0) throw("no circle thru colinear points");
    det = 0.5/det;
    asq = a0*a0 + a1*a1;
    csq = c0*c0 + c1*c1;
    ctr0 = det*(asq*c1 - csq*a1);
    ctr1 = det*(csq*a0 - asq*c0);
    rad2 = ctr0*ctr0 + ctr1*ctr1;
    return Circle(Point<2>(ctr0 + b.x[0], ctr1 + b.x[1]), sqrt(rad2));
}
// triangular element
struct Triel {
    Point<2> *pts;
    int p[3];
    int d[3];
    int stat;
    void setme(int a, int b, int c, Point<2>* ptss)
    {
        pts = ptss;
        p[0] = a; p[1] = b; p[2] = c;
        d[0] = d[1] = d[2] = -1;
        stat = 1;
    }
    int contains(Point<2> point)
    {
        double d;
        int i, j, ztest = 0;
        for (i = 0; i < 3; i++) {
            j = (i+1)%3;
            d = (pts[p[j]].x[0] - pts[p[i]].x[0])*(point.x[1]-pts[p[i]].x[1]) -
                (pts[p[j]].x[1] - pts[p[i]].x[1])*(point.x[0]-pts[p[i]].x[0]); 
            if (d<0.0) return -1;
            if (d==0.0) ztest= 1;
        }
        return (ztest? 0:1);
    }
};

double incircle(Point<2> d, Point<2> a, Point<2> b, Point<2> c)
{
    Circle cc = circumcircle(a, b, c);
    double radd = sqrt(d.x[0] - cc.center.x[0]) + sqrt(d.x[1]-cc.center.x[1]);
    return (sqrt(cc.radius) - radd);
}


struct Delaunay {
    int npts, ntri, ntree, ntreemax, opt;
    double delx, dely;
    vector<Point<2>> pts;
    vector<Triel> thelist;
    vector<list<int>> belong;

    int* perm; //permutation
    Delaunay(vector<Point<2>> & pvec);
    void insertapoint(int r);
    int storetriangle(int a, int b, int c);
    int whichcontainspt(const Point<2> &p);
    void erasetriangle(int a, int b, int c, int d0, int d1, int d2);
};

Delaunay::Delaunay(vector<Point<2>> & pvec)
{
    npts = pvec.size();
    ntri = 0;
    ntree = 0;
    pts.reserve(npts+3);
    thelist.reserve(2*npts);
    for (int j = 0; j<npts+3; j++) {
        list<int> lst = {0};
        belong.push_back(lst);
    }
    ntreemax = 10*npts + 1000;
    double xl, xh, yl, yh;
    perm = new int[npts];
    xl = xh = pvec[0].x[0]; 
    yl = yh = pvec[0].x[1]; 
   
    for (int j = 0; j<npts; j++) {
        pts[j] = pvec[j];
        perm[j] = j;
        if (pvec[j].x[0] < xl) xl = pvec[j].x[0];
        if (pvec[j].x[0] > xh) xh = pvec[j].x[0];
        if (pvec[j].x[1] < yl) yl = pvec[j].x[1];
        if (pvec[j].x[1] > yh) yh = pvec[j].x[1];
    }
    
    delx = xh - xl;
    dely = yh - yl;
    pts[npts] = Point<2>(0.5*(xl+xh), yh+dely);
    pts[npts+1] = Point<2>(xl-0.5*delx, yl-0.5*dely);
    pts[npts+2] = Point<2>(xh+0.5*delx, yl-0.5*dely);
    storetriangle(npts, npts+1, npts+2);

    for(int j = 0; j < npts; j++) insertapoint(perm[j]);
}

void Delaunay::insertapoint(int r)
{
    int i,j,k,l,s,tno,ntask, d0, d1, d2;
    unsigned long key;
    int tasks[50];
    int taski[50];
    int taskj[50];
    for(j = 0; j < 3; j++) {
        tno = whichcontainspt(pts[r]);
        if (tno >= 0) break;
    }
    i = thelist[tno].p[0];
    j = thelist[tno].p[1];
    k = thelist[tno].p[2];

    d0 = storetriangle(r,i,j);
    tasks[++ntask] = r; taski[ntask] = i; taskj[ntask] = j;
    d1 = storetriangle(r,j,k);
    tasks[++ntask] = r; taski[ntask] = j; taskj[ntask] = k;
    d2 = storetriangle(r,k,i);
    tasks[++ntask] = r; taski[ntask] = k; taskj[ntask] = i;
    erasetriangle(i,j,k,d0,d1,d2);

    while(ntask) {
        s = tasks[ntask]; i = taski[ntask]; j = taskj[ntask--];
        // find edge corrspond point l
        for(auto it = belong[i].begin(); it != belong[i].end(); it++)
        {
            int l = 0;
            Triel tt = thelist[*it];
            for(int q = 0; q < 3; q++)
            {
                if(tt.p[q] == i)
                {
                    if(tt.p[(q+1)%3] != s && tt.p[(q-1)%3] == j)
                    {
                        l = tt.p[(q+1)%3];
                        break;
                    }
                    else if (tt.p[(q-1)%3] != s && tt.p[(q+1)%3] == j)
                    {
                        l = tt.p[(q-1)%3];
                        break;
                    }
                }
            }
        }
        if(incircle(pts[l], pts[j], pts[s], pts[i]) > 0.0)
        {
            d0 = storetriangle(s,l,j);
            d1 = storetriangle(s,i,l);
            erasetriangle(s,i,j,d0,d1,-1);
            erasetriangle(l,j,i,d0,d1,-1);
            tasks[++ntask] = s; taski[ntask] = l; taskj[ntask] = j;
            tasks[++ntask] = s; taski[ntask] = i; taskj[ntask] = l;
        }
    }
}
int Delaunay::whichcontainspt(const Point<2> &p)
{
    int i,j,k = 0;
    while (thelist[k].stat <= 0)
    {
        for (i = 0; i < 3; i++)
        {
            if ((j = thelist[k].d[i]) < 0)
                continue;
            if (thelist[j].contains(p) > 0)
                break;
        }
        if (i == 3)
            return -1;
        k = j;
    }
    return k;
}

int Delaunay::storetriangle(int a, int b, int c)
{
    thelist[ntree].setme(a, b, c, &pts[0]);
    belong[a].push_back(ntree);
    belong[b].push_back(ntree);
    belong[c].push_back(ntree);
    if (++ntree == ntreemax)
        throw("thelist is too small");
    ntri++;
    return ntree-1;
}

void Delaunay::erasetriangle(int a, int b, int c, int d0, int d1, int d2)
{
    int j;
    for (auto it = belong[a].begin(); it != belong[a].end(); ++it) {
        int match = 0;
        for(int q = 0; q < 3; q++)
        {
            if(thelist[*it].p[q] == b)
                match++;
            if(thelist[*it].p[q] == c)
                match++;
        }
        if (match == 2)
        j = *it;
    }
    unsigned long key;
    
    thelist[j].d[0] = d0;
    thelist[j].d[1] = d1;
    thelist[j].d[2] = d2;
    thelist[j].stat = 0;
    ntri--;
}



int main(void)
{
    Point<2> a = {0.0, 1.0};
    Point<2> b = {1.0, 1.0};
    Point<2> c = {1.0, 0.0};

    vector<Point<2>> d;
    d.push_back(a);
    d.push_back(b);
    d.push_back(c);
    Delaunay asd = Delaunay(d);
}